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Abstract. An accurate non-gradient-expansion based correction to Thomas-Fermi 
is developed using solvable model. The used model is a system of N non-interacting 
electrons moving independently in the Coulomb field of the nuclear charge. The 
presented correction is applicable for atoms and should be extendable beyond that. The 
method exploits the fact that the difference between the Thomas-Fermi approximation 
and the non-interacting kinetic energy is comparable to the difference between the 
same values inside the proposed solvable model. The numerical experiments show 
that by adding this correction factor, the precision of Thomas-Fermi approximation is 
enhanced by an order of magnitude. 


1. Introduction 

As known, the original density fnnctional theory (DFT) is based on Hohenberg and 
Kohn prominent work [T] where they proved that the ground state of any many-electron 
system is completely characterized by its density and that the energy functional of the 
system attains its minimum at the density corresponding to the ground state. However, 
representing the contribution of kinetic energy as density functional has proven to be 
challenging as the accuracy and applicability of the proposed kinetic energy density 
functionals (KEDF) are generally not sufficient [2], [3l HI |5]. Alternatively, Kohn and 
Sham (KS-DFT) [6] suggested an approach where the “orbitals” are reintroduced where 
the sum of the orbitals’ densities equals to the density of the real system and the 
kinetic energy is defined as the kinetic energy of the introduced “fictitious” system. 
Computational-wise, this results in converting the problem back from 3-dimensional 
(3D) to 3N dimension as N orbitals are determined by solving the governing N equations 
self-consistently mmn, where N is the number of the particles. Despite this drawback, 
KS-DFT nowadays dominates atomistic calculations [7j. In parallel, the work to find 
an accurate KEDF is still active, but modestly. In the recent years, it starts gaining 
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attention To distinguish this DFT doctrine from KS-DFT which depends 

on the orbital, it is commonly named “orbital-free’ DFT (OF-DFT) |T2l[T3]. 

Since the first independently proposed KEDF by Thomas [ll] and Fermi [15] (TF), 
a huge number of KEDFs has been suggested. However, Thomas-Fermi model with 
various corrections dominates the field [21 El IS] • The usual corrections are either based 
on gradient expansion [T6l ITT] [T8] or on adding additional class of KEDF P HH [20l |2T] . 
Furthermore, Thomas-Fermi-based KEDF are used in some applications satisfying the 
assumption that the density is nearly uniform. For example, they are used for metals 123 
and warm dense matter [231 [22] . As for the correction, it would be assumed abstractly 
that gradient expansion should pave a reasonable route. However, it is known that the 
high order gradient-based corrections diverges for finite systems. Therefore, seeking non¬ 
gradient-expansion based correction is desirable. One of the recent corrections suggested 
by Burke and coworkers [5] is based on uniform WKB analysis in 1-dimensional case. 
The addition has neither sums nor derivatives. 

In this paper, we present a new non-gradient based method to enhance the Thomas- 
Fermi approximation by adding a correction factor derived using an exactly solvable 
model. The used model is for N non-interacting particles moving in the Coulombie 
field of the nuclear charge. This allows using the simple Rydberg formula to calculate 
part of the energy and the particle density is expressed in an analytic form through 
Laguerre polynomials and exponential functions. The correction exploits the fact that 
the difference between the Thomas-Fermi approximation and the non-interacting kinetic 
energy is comparable to the difference between the same values inside the proposed 
solvable model. The concept can be presented best using the following equation 

Ts = Ttf + (Ts — Ttf) ~ Ttf -I- (T — Ttf)- (1) 

where Tg is the non-interacting kinetic energy and T^f is Thomas-Fermi approximation. 
In the same equation we use the notation T and Ttf for the kinetic energy and 
Thomas-Fermi energy of the proposed model. Our work follows a similar approach that 
was developed to estimate the correlation energy for two-electron atoms [25]. In our 
numerical experiments, we show that the correction factor can increase the precision of 
Thomas-Fermi approximation by around an order of magnitude for atoms. To apply it 
for molecules, the method needs further extensions. As shown in Q, the kinetic energy 
is expressed as Thomas-Fermi energy plus the correction obtained from the exactly 
solvable model. We calculated the results numerically for a wide range of atoms up to 
xenon (Z = 54). Our results appear to be more accurate in comparison with Thomas- 
Fermi approximation by a factor between 11 (for helium) and 72 for xenon. We also 
discuss the large-Z limit of the energy and the density and compare the correction to 
those by gradient expansions. 

2. The correction 

As aforementioned, the general idea of the proposed correction is that the difference 
between the Thomas-Fermi approximation and the non-interacting kinetic energy is 
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comparable to the difference between the same valnes inside the proposed solvable model 
as shown in Q. This follows a snccessfnl and analogons method used to estimate the 
correlation energy for two-electron atoms iza . Although we assume the model of an 
N non-interacting particles moving in the Coulombic field of the nuclear charge, the 
method can be applied for a wide range of potential models. In this paper, the focus is 
to present the approach and to illustrate its applicability for simple systems like atoms. 
So, we suggest the use of a model that in many respects resemble the atoms and that 
has an analytical quantum-mechanical solution. Another important reasons for selecting 
such a model is that it is possible to write the correction T — Ttf as a function of only 
the electrical charge Z. Furthermore, this form for T and Ttf provides more in depth 
understanding of the the model as it will be shown. 

Since we are interested in an equivalent system of non-interacting particles, the 
energy T does not include correlation effects. The inter-electron interaction affects 
the effective potential and therefore the only difference between our model and an 
atom is disregarding the screening of nuclear charge by inner electrons. The presented 
model retains Coulomb singularity at the origin, as well as Coulomb attraction for large 
distance, however the attraction force at large distance is much larger in the presented 
model. The region of applicability of the Thomas-Fermi model, dp~^/dr ^ 1, where 
Pfj_ = y^2(/i — V) and p is the chemical potential, is violated in a small region of radius 
To = Z~^ adjoining the nucleus, where quantum effects become significant. In the 
neighborhood of the nucleus, the field is practically identical to the Coulomb field —Z/r. 
In our approach however, we treat the Coulomb problem exactly, without considering 
the quantum corrections separately. 

In the assumed model, we consider a system with the nuclear charge Z = N. 
Thus, the screened Coulomb potential for the atom has the same behavior —Z/r at 
the origin, but differs far from it. Without inter-particle interaction, each electron can 
be considered as occupying an orbital characterized by the principal quantum number 
n = 1,2,.... For example, two electrons with n = 1 in a configuration form the 
closest inner K-shell (in X-ray notations), eight electrons with n = 2 form L-shell in a 
configuration 2s^2p®. Generally, a completely filled n-th shell has 2n‘^ electrons with 
possible quantum numbers Z = 0,1,..., n — 1, m = —I, + and spins a = ±|. 

Let us denote Umax as the last shell with a non-zero occupation number. In the 
ground state, all lowest shells with n < Umax are completely filled, while the last shell 
with n = Umax can be filled either partially or completely. For simplicity, we consider 
initially only the states in closed-shell configurations, i.e. when all shells up to n = Umax 
are completely filled. The same model can be extended for other configurations either 
by considering angular dependence or simply by interpolation between closed shells as 
it will be shown later. The total number of electrons on shells with n = 1,2,, Umax 
can be obtained by summation of the occupation numbers for each individual shell, 

^max 

N = ^ 2n^ = -nmax(nmax + l)(2nmax + !)• (2) 

n=l 

The kinetic energy of the proposed model can be calculated by exploiting the kinetic 
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energies for individual orbitals. This energy for each orbital is given by Rydberg formula, 
that is Z^/(2n^) in atomic units, where n is the principal quantum number. Now, the 
total kinetic energy can be calculated by summation of contributions from each shell, 

T = Z 2"' (3) 

where the kinetic energy of the model T corresponds to the same value in (^. 

To fully specify the proposed functional, we also need to calculate the Thomas- 
Fermi approximation of energy for the model. As previously stated, the Ttf can be 
presented as a function of the electrical charge Z when the proposed model is used. 
For closed-shell conhguration which are spherically symmetric, Thomas-Fermi kinetic 
energy is 

poo 

Ttf = 4:71 rVo(r)dr, (4) 

Jr=0 

where the KEDF tq is 

Mr) = ^ (37r^y^%^/^{r). (5) 

To be able to calculate the Thomas-Fermi approximation using Q, for the proposed 
model, it is necessary to have the corresponding electron density p. From the definition 
of the model, we know that p can be acquired by combining the wave functions for 
individual electrons, as follows. A wavefunction of an electron on the orbital (n, Z, m) is 

= Rn,l{r)^l,m{0,(t)), ( 6 ) 


where is a wave function, n is the principal quantum number, I is an orbital quantum 
number, m is an azimuthal quantum number, r is the radius, Yi^rn is a spherical harmonic, 
6^ is a polar angle, 0 is an azimuthal angle. In (|^, Rn^i represents the the radial 
component of the wavefunction and it is defined using the following equation 


Rn.iir) = , — 


(2Z\^ (n-Z-1)! 


Zr\ /2ZrV 


/2Zr\ 


\ [-) Mirror “p [--) (—) (—j ■ (7) 


Using ([^ the electronic density of an atom p in our model can be given 

^max CL 1 I 

m = E E E 


as 


( 8 ) 


n=l Z=0 m=—l 

By exploiting the properties of the spherical harmonics and by summing over m, the 
above equation is reduced to 

^max CL 1 

(9) 


1 


pW = z;EE(2' + iXiW- 

n=l 1=0 


which is angular independent. Finally we can calculate the value of Ttf by incorporating 
p, given by equations ([^ and (^. 
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After determining the values for T and we can explicitly write the functional 
for calculating the kinetic energy for an atom using the correction based on the proposed 
model. The new functional is given by the following equation 

T ^ Ttf + 5T, (10) 

where dT is found by considering the exactly solvable model with the same nuclear 
charge Z, 

6f = f- Ttf- (11) 

When using the proposed functional for the kinetic energy of an atom, we would 
calculate the Ttf based on some electron density approximation and the correction 
using the value of Z. In practical application, it is necessary to consider the case when 
the last shell is partially hlled. In such cases, it is not possible to directly apply (|^ with 
an integer Umax- Since the closed shells occur only for a few values of the nuclear charge 
given by the sequence of “magic numbers” 2,10, 28, 60,110,..., we need to dehne the 
interpolation of the function Z i—)■ dT' to other integer values of Z. Here, we use an 
interpolation by a cubic polynomial based on four points Z = 2, 10, 28, and 60: 

df = 0.21210 - 0.19860Z + 0.12815^^ + O.OOOlOZl (12) 


3. The Z expansion of the model 


For many applications, it is more convenient to represent the energy of the system based 
on the total nuclear charge Z. This form can be acquired by solving (|^ in respect to 
U max and substitution of the result into ([^, we obtain the energy as a function of the 
nuclear charge. 


T = 


^ 3 _i/ 3 ^_i ^ 3 - 2 / 3 ^ - 1 ) D= {bAZ + V2916Z2 - 3 ) .(13) 

A more suitable format for (13) is in the form of an expansion in powers of Z~^l‘^. The 
transformed equation has the following form 


f ~ (3/2) 


H_ ^ v5/3_ \ _^1/3 

2 6 X 12V3 3888 x 18^3 


1 

^ 69 984 X 12V3 


^-1/3 + C)(^-5/3)^ 


(14) 


where the terms proportional to and Z^ are identically zero. In 

numerics, (14) is simply 


T ~ 1.144714Z^/^ - 0.5^2 + 0.072798Z^/^ 

- 0.000098Z^/3 ^ o.000006Z-^/3 + 0{Z-^l^). (15) 


With the goal of being able to assess the validity of the model for calculating 
the kinetic energy, we compare it to the corresponding expansion of Thomas-Fermi 
approximation, where the energy of an atom with large nuclear charge Z can be written 








Correction to kinetic energy density using exactly solvable model 


6 


as an asymptotic series in powers of a small parameter Z This expansion is found 
to be 


T ~ 0.7687452'''^ - 0.52'^ + 0.2699002’'''“, 


(16) 


The largest absolute difference between equations (14) and (16) is the leading term. This 


is mainly due to the inclusion of Coulomb repulsion between electrons which decreases 
the Thomas-Fermi energy by a factor of 1.489 because of the increase of the size of the 
atom. The second subdominant term ~ Z‘^ remains the same. This term comprises the 
correction of strongly bound electrons [271 EH] which is not affected by the inter-electron 
repulsion. 

To assess the accuracy of the Thomas-Fermi approximation, we analyze the 
asymptotic behavior of the Thomas-Fermi energy in the limit of large number of 
electrons N = Z. To estimate the large-A^ behavior, we have calculated Ttf for 
increasing values of Umax and found, using Richardson’s extrapolation [29], that for 
large Umax 


Tr 


TF 


1.1447142“''“ - 0.6258562“ + 0.1468782“''“. 


(17) 


We characterize the accuracy of the approximation by measuring the deviation of the 
expansion ( pTj ) from the exact coefficients in (16) [30] . We find that the Thomas-Fermi 
approximation correctly reproduces the leading term 1.144714Z^/^, but makes a 25% 
error in the subdominant term —0.5Z^. Thus we can say according to the definition 
from [31] that the Thomas-Fermi approximation is large-iV asymptotically exact to the 
first degree (AEl). 


4. Numerical results and discussion 

In this section, the presented functional and the corresponding approach are 
implemented for verification. First we compare the electron density of the proposed 
model and its large-Z limit. Then, we compare the model to the known gradient 
expansion corrections of Thomas-Fermi KEDF. Finally, we perform computational 
experiments to show that the proposed functional gives a significant improvement to 
the standard Thomas-Fermi approximation when applied to atoms. An improvement 
by about an order of magnitude is achieved. 

4-1. Electron density in the limit of large Z 

In this subsection, we analyze the behavior of electron density of the proposed model 
to the one acquired using the Thomas-Fermi formalism vs. Z. We start with some 
remarks regarding the density of the model. In Section the electron density is derived 
from the wavefunction as a sum of nmax(^^max + l)/2 terms given by 0. For large Z, 
the number of terms grows to infinity. Here, we use an alternative approach based on 
Thomas-Fermi formalism to derive the limit of the density for large values of Z in a 
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more explicit way. Without an inter-particle interaction, the Thomas-Fermi equation 
relating the electron density and the potential takes especially simple form|32j. 

1 


where fi can be determined from the equation (assuming spherical symmetry) 
/ Anr p{r) dr = N. 


(18) 


(19) 


In (19), N is the number of electrons and is the Thomas-Fermi radius of the atom. 


or a turning point in the potential 17, that is determined from the equation 

V{r^)=p. (20) 

For the Coulomb potential V = —Z/r, we obtain 

Z , . 

rm =- 21 

F 


From equations (19) and (21), it can be shown that 

r^ = Z-^{?,V2N)''\ ( 22 ) 

As we consider the the case of a neutral atom [N = Z), then, the dependence on 
Z can be eliminated by introducing scaled radius and scaled electron density 

f = Z'/’r, (23) 

P = P/Z\ (24) 


SO that now (18) can be rewritten as 


2^/2 n 


-1/3 


3/2 


(25) 


(26) 


if f < fm and zero otherwise, where the scaled turning point is 
rm = 18'/^ 

Figure [^shows the scaled densities p{r) = Z~‘^p{Z^^^r), where p is given by ([^ for 
increasing numbers of electronic shells, rimax = 1,2, 3, 5, calculated using ([8| together 


with the limiting case given by (25). It is clear that the deviation from the Thomas- 
Fermi limit, p — ptf is an oscillating function having exactly Umax local maxima, see 
the insert in Figure The amplitude of the oscillations decrease as Z increases. 
These oscillation effects are related to the shell structure, with maxima corresponding 
to the filled shells. Clearly, the Thomas-Fermi model describes only the averaged 
physical quantities, and so it requires a special generalization [33] to treat such spatial 
irregularities. In the recent paper of Burke and coworkers [5], their approximation 
improves the accuracy everywhere including the turning points and reproduces the 
oscillations. 
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^max ~ ^ 
^max ~ 2 
^max ~ 3 
^max ~ 5 


0.0 0.5 1.0 1.5 2.0 2.5 3.0 

r 

Figure 1. Scaled densities for increasing values of Z and Thomas-Fermi limit at Z ^ oo. 



4-2. Thomas-Fermi vs. gradient expansion for the exactly solvable model 


The leading term of gradient expansion for the kinetic energy is just the Thomas-Fermi 
energy given by equations (|^ and We use the notation Tq = = Ttf- The 

second order term is dehned through Weizsacker correction Tw as 
1 (o' 

T2 = -Tw, Tw = dvr / dr, ( 27 ) 

9 Jr=0 op 

and the second order approximation is = Tq + T 2 , where spherical symmetry is 
assumed. In fourth order of the gradient expansion, we have 

POO 

r4 = dvr / r^ r4(r) dr, = Tq + T2 + T4, ( 28 ) 

Jr=0 

where 


r4(f) 


(3-^u/ 3 
5d0 ^ 



9(2^ + p")(p')2 1(p')4 

8 ■ ' ' ^ 3 


(29) 


There were several studies of the gradient expansion for atoms [Ml l35| IM] , The 
Thomas-Fermi method always underestimates the energy, and the accuracy slowly 
improves with increase of number of electrons, remaining on the level of few percent 
even for heavy atoms. The hrst correction of gradient expansion always improves the 
accuracy, but applying the fourth order correction of gradient expansion typically makes 
the results worse. The analysis of trends of the gradient expansion for atoms remains 
somehow inconclusive, because of very slow asymptotical behavior (typically as -1/3) 
and because of shell effects (oscillations of density). 

For the exactly solvable model, relative error as a function of the number of shells 
is shown on Figure where the non-interacting kinetic energy Tg = T is dehned by 
(|^, the Thomas-Fermi energy is dehned by (|^, and the gradient expansion is dehned 
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Figure 2. Relative error of Tg for Thomas-Fermi = Ttf, marked by circles), second 

(r(2) = To + T 2 , square markers), and fourth order approximation (T^^^ = Tq + T 2 + T 4 , 
diamond markers) in gradient expansion. The relative errors are defined as (Tg — 
where n = 0 , 2 ,4. 


by equations ( |27[) and (28). For the exactly solvable model, the relative error to Tg is 
shown on Figurej^as a function of the number of shells. The values for are 

calculated by substituting p into the corresponding equations. For large Z, including the 
second order correction clearly improves the accuracy by a factor of 6, and including the 
fourth order correction improves the accuracy by an additional factor of 3. However, this 
trend is visible only for large Umax, starting from Umax ~ 10 corresponding to unrealistic 
Z > 1000. It explains the fact that the forth order correction for atoms improves the 
accuracy only for very heavy atoms. 

Asymptotic behavior of accuracy at large Z is shown on Figure]^ The asymptotical 
expansion for the Thomas-Fermi energy has been presented in the previous section and 
has the form given in dl?] ). In a similar way the terms of gradient expansion were 
determined numerically and have the following forms 


T 2 ~ 0.10942^2 + 0.045Z®/3 + .... 


Notice that the expansion in (30) starts from the term Z^ 


i.e. 


(30) 

the coefficient of the 


leading term ~ is zero. This statement can be proven rigorously too, by calculating 
the integral (27) using the same semiclassical approximation. A similar expansion 
(starting from the term ~ Z^) was found for the fourth-order correction to the kinetic 
energy, 


Ti ~ 0.015052Z2 + 0.0078Z'^/^ + .... (31) 

These calculation show that the leading term of the exact energy, 1.447Z^/^, is 
absorbed by the leading term of Thomas-Fermi energy, and the subdominant term of 
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Figure 3. Relative error in logarithmic scale for very large number of electronic shells for 
Thomas-Fermi approximation and the approximations including second and fourth order 
terms of gradient expansion. The markers are the same as on Figure Since the asymptotic 
dependence is linear with the tangent (-1) in all cases, it shows that the relative errors at 
large Z are proportional to in all cases. 


the exact energy, ~ —0.5Z^, can be accurately resumed from the corresponding terms 
in the series = Ttf + T 2 + T 4 . 


4^.3. Calculations for atoms 


In this subsection we discuss the effectiveness of the functional given by (10) for 


calculating the kinetic energy of atoms. This is done by calculating the kinetic energy 


for atoms using (10) and compare the results with Hartree-Fock energies. To show 


explicitly the advantages of the proposed functional we compare the kinetic energies 
calculated using the new approach with the Thomas-Fermi approximation and the 
standard functionals based on gradient expansion. We calculated Thomas-Fermi kinetic 
energy for atoms with closed shells using the electronic densities derived from dementi 


Tables [37]. The correction 6T was calculated by (11) in case of closed shell atoms and 
by (12) in partially hlled cases. The results are shown in Table 

We found that including the correction 6T increases the accuracy of Thomas-Fermi 
approximation by more than nine times. As an approximation to the energy Tg, here we 
used the Hartree-Fock energies tabulated in [37] . The results are especially accurate for 
atoms with spin-paired electrons (He, Be, Ne, Mg, Ar, Kr, and Xe), where our numerical 
results always give the upper bound for the kinetic energy, while the Thomas-Fermi 
approximation gives lower bounds in all cases. We expect that spin-polarized version 
of Thomas-Fermi theory would work signihcantly better for cases with uncompensated 
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Table 1. Comparison of Thomas-Fermi kinetic energy with the improved Thomas-Fermi, 
for atoms. Results from the second and fourth order gradient expansion are included too. 
The relative error is calculated in comparison with Hartree-Fock kinetic energies taken from 
dementi Tables m which are known to be close to the energy Tg. The density for atoms 
with nonzero angular momentum is spherically averaged. 


z 

Atom 

Ttf 

Relative error, % 

Ttf + T 2 Ttf + T2 -|- T4 

Ttf + Sr 

2 

He 

-11 

0.59 

3.6 

0.95 

3 

Li 

-10 

0.62 

3.1 

0.26 

4 

Be 

-9.9 

0.50 

2.9 

0.21 

5 

B 

-10 

-0.53 

1.6 

-0.52 

6 

C 

-11 

-1.3 

0.67 

-1.0 

7 

N 

-11 

-1.7 

0.16 

-1.2 

8 

0 

-10 

-1.6 

0.10 

-0.95 

9 

F 

-9.4 

-1.2 

0.37 

-0.47 

10 

Ne 

-8.4 

-0.56 

0.95 

0.28 

11 

Na 

-8.1 

-0.49 

0.94 

0.37 

12 

Mg 

-7.8 

-0.44 

0.95 

0.43 

13 

AI 

-7.6 

-0.45 

0.89 

0.42 

14 

Si 

-7.5 

-0.47 

0.83 

0.39 

15 

P 

-7.4 

-0.50 

0.77 

0.36 

16 

s 

-7.3 

-0.51 

0.72 

0.35 

17 

Cl 

-7.1 

-0.51 

0.70 

0.35 

18 

Ar 

-7.0 

-0.49 

0.69 

0.36 

36 

Kr 

-5.8 

-0.69 

0.18 

0.11 

54 

Xe 

-5.2 

-0.68 

0.067 

0.073 


spins, like atoms of B, C, N, O, F. 

Overall, the same table shows that the proposed method is competitive, with 
functionals based on the second and forth order gradient expansions when precision is 
considered. It manages to hnd approximations with a lower error than the second/forth 
order expansion for 18/14 out of 19 tested atoms. From the tabulated results, it is 
noticeable that the new functional is more robust than the other functionals in the 
sense that it has a good performs for both small and large atoms. It is important 
that the proposed correction procedure is much simpler to calculate, since it is a single 
variable function compared to very complicated electron density functionals in case of 
the gradient based corrections. 
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In this paper we have presented a non-gradient-based correction to Thomas-Fermi 
fnnctional for atoms. The presented approach is general and should be extendible for 
molecules. The method uses an auxiliary system of noninteracting electrons that is in 
many respects similar to the atomic system with the same number of electrons. It results 
in simplifying the calculations considerably when compared to the one based on gradient 
expansion. The obtained accuracy is improved by more than nine times in comparison 
to Thomas-Fermi model. Our numerical test have also shown that the proposed method 
manages to achieve similar, slightly better, precision that the standard gradient based 
functionals. As for the density, our presented approach allows the characteristic shell 
oscillations. 

This type of approach can potentially be extended to systems other than atoms. 
One example is the modeling on A^-electron quantum dots, where we could consider 
another solvable model, of N non-interacting particles bound in a harmonic potential. 
Our initial test on this problem have given promising results. 
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